// copyright:no copyright,use it if you like
// license boilerplate:BSD 3-Clause
// filename:Interpolation.hpp
// author:Tempest(https://gitee.com/Tempest56890)
// date:2022-10-20
// brief:implement the Newton formula and Hermite fomula.

#ifndef __TEMPEST_INTERPOLATION__
#define __TEMPEST_INTERPOLATION__

#include <iostream>
#include <cmath>
#include <limits>
#include <vector>
#include <stdlib.h>

const double _eps=std::numeric_limits<double>::epsilon();

// implement the funciton of BaseClass.
// improve a virtual function operator().
// example:
//    (Temporarily omitted.)
class Function {
public:
    // Defines a virtual function operator().
    virtual double operator()(const double x) const =0;
    // First-Order derivative function is defined.
    virtual double First_Order_diff(const double &x) const {
        return ((*this)(x+100*_eps)-(*this)(x-100*_eps))/(_eps*200);
    }
    // Second-Order derivative function is defined.
    virtual double Second_Order_diff(const double &x) const {
        double left_diff=((*this)(x)-(*this)(x-200*_eps))/(_eps*200);
        double right_diff=((*this)(x+200*_eps)-(*this)(x))/(_eps*200);
        return (right_diff-left_diff)/(_eps*200);
    }
};

// implement the funciton of BaseClass.
// implements derivative functions.
// example:
//    (Temporarily omitted.)
class Polynomial {
protected:
    int order;
    // coef is used to store the coefficients of the polynomial。
    std::vector<double> coef;
public:
    Polynomial() : coef({0.0}), order(0) {};
    Polynomial(std::vector<double> coef,int order)
    :coef(coef),order((int)coef.size()-1) {}
    // generate a simple polynomial.
    Polynomial(const double x0)
    {
        order = 0;
        coef.resize(1);
        coef[0]=x0;
    }
    // generate a simple polynomial.
    Polynomial(const double x0, const double x1)
    {
        order = 1;
        coef.resize(2);
        coef[0]=x0;
        coef[1]=x1;
    }
    double operator() (const double x) const {
        double sigma=0.0,xx=1;
        for(int i=0;i<=order;i++) {
            sigma+=xx*coef[i];
            xx*=x;
        }
        return  sigma;
    }
    // refactoring + Operator
    Polynomial operator+ (const double rhs) {
        Polynomial res=*this;
        res.coef[0]+=rhs;
        return res;
    }
    // refactoring + Operator
    Polynomial operator+ (const Polynomial& rhs) {
        Polynomial res=*this;
        int _order=(order>=rhs.order ? order : rhs.order);
        if(order>0&&order==rhs.order&&coef[order]+rhs.coef[order]==0.0) {
            _order-=1;
        }
        res.coef.resize(_order+1);
        for(int i=0;i<=std::min(order,rhs.order);i++) {
            res.coef[i]+=rhs.coef[i];
        }
        if(order<rhs.order){
            for(int i=order+1;i<=rhs.order;i++) {
                res.coef[i]=rhs.coef[i];
            }
        }
        res.order=(int)res.coef.size()-1;
        return res;
    }
    // refactoring - Operator
    Polynomial operator- (const double rhs) {
        Polynomial res=*this;
        res.coef[0]-=rhs;
        return res;
    }
    // refactoring - Operator
    Polynomial operator- (const Polynomial& rhs) {
        Polynomial res=*this;
        int _order=(order>=rhs.order ? order : rhs.order);
        if(order>0&&order==rhs.order&&coef[order]-rhs.coef[order]==0.0) {
            _order-=1;
        }
        res.coef.resize(_order+1);
        for(int i=0;i<=std::min(order,rhs.order);i++) {
            res.coef[i]-=rhs.coef[i];
        }
        if(order<rhs.order){
            for(int i=order+1;i<=rhs.order;i++) {
                res.coef[i]=-rhs.coef[i];
            }
        }
        res.order=(int)res.coef.size()-1;
        return res;
    }
    // refactoring * Operator
    Polynomial operator* (const double rhs) {
        Polynomial res;
        if(rhs!=0.0) {
            res.coef.resize(order+1);
            for(int i=0;i<=order;i++) {
                res.coef[i]=coef[i]*rhs;
            }
        }
        res.order=(int)res.coef.size()-1;
        return res;
    }
    // refactoring * Operator
    Polynomial operator* (const Polynomial& rhs) {
        Polynomial res;
        res.order=order+rhs.order;
        res.coef.resize(res.order+1);
        for(int i=0;i<=order;i++){
            for(int j=0;j<=rhs.order;j++){
                res.coef[i+j]+=coef[i]*rhs.coef[j];
            }
        }
        res.order=(int)res.coef.size()-1;
        return res;
    }
    // calculate the derivative of a polynomial and output it in Polynomial class format。
    Polynomial diff() const {
        Polynomial res=*this;
        if(order==0) {
            res.coef[0]=0.0;
        } else {
            res.coef.resize(order);
            for(int i=0;i<=order-1;i++) {
                res.coef[i]=coef[i+1]*(i+1);
            }
        }
        res.order=(int)res.coef.size()-1;
        return res;
    }
    // refactoring << Operator
    friend std::ostream& operator<<(std::ostream& out, const Polynomial& pol) {
        if(pol.coef[pol.order]!=0.0) {
            out<<pol.coef[0];
        }
        for(int i=1;i<pol.order;i++) {
            if(pol.coef[i]<0&&fabs(pol.coef[i])!=1.0) {
                out<<pol.coef[i]<<"*x.^"<<i;
            } else if (pol.coef[i]>0&&fabs(pol.coef[i])!=1.0) {
                out<<"+"<<pol.coef[i]<<"*x.^"<<i;
            }
        }
        if(pol.order!=0) {
            if(pol.coef[pol.order]<0&&fabs(pol.coef[pol.order])!=1.0) {
                out<<pol.coef[pol.order]<<"*x.^"<<pol.order<<std::endl;
            } else if (pol.coef[pol.order]>0&&fabs(pol.coef[pol.order])!=1.0) {
                out<<"+"<<pol.coef[pol.order]<<"*x.^"<<pol.order<<std::endl;
            }
        }
        return out;
    }
};

// implement the funciton of Derived class.
// implements Newton Interpolation.
// example:
//    (Temporarily omitted.)
class Newton {
private:
    std::vector<double> x,y;
    // used to store the difference table.
    std::vector<std::vector<double>> diff_table;
public:
    Newton(const std::vector<double>& x,const std::vector<double>& y)
    :x(x),y(y) {
        diff_table.resize(x.size());
        for(int k=0;k<=x.size()-1;k++) {
            diff_table[k].resize(x.size());
            diff_table[k][0]=y[k];
        }
        for(int i=1;i<x.size();i++) {
            for(int j=1;j<=i;j++) {
                diff_table[i][j]=(diff_table[i][j-1]-diff_table[i-1][j-1])/(x[i]-x[i-j]);
            }
        }
    }
    // generate a polynomial by Newton Interpolation.
    Polynomial getPol() const {
        Polynomial res,tem(1);
        for(int i=0;i<=x.size()-1;i++){
            res=res+tem*diff_table[i][i];
            tem=tem*Polynomial(-x[i],1);
        }
        return res;
    }
};

// implement the funciton of Derived class.
// implements Hermite Interpolation.
// example:
//    (Temporarily omitted.)
class Hermite {
private:
    //d_o represents the derivative order of y corresponding to x.
    std::vector<double> x,y,d_o;
    std::vector<std::vector<double>> diff_table;
public:
    Hermite(const std::vector<double>& x,const std::vector<double>& y,const std::vector<double>& d_o)
    :x(x),y(y),d_o(d_o) {
        diff_table.resize(x.size());
        for(int k=0;k<=x.size()-1;k++) {
            diff_table[k].resize(x.size());
            diff_table[k][0]=y[k-d_o[k]];
        }
        for(int i=1;i<x.size();i++) {
            for(int j=1;j<=i;j++) {
                if(x[i]!=x[i-j]) {
                    diff_table[i][j]=(diff_table[i][j-1]-diff_table[i-1][j-1])/(x[i]-x[i-j]);
                } else {
                    diff_table[i][j]=y[i];
                }
            }
        }
    }
    // generate a polynomial by Hermite Interpolation.
    Polynomial getPol() const {
        Polynomial res,tem(1);
        for(int i=0;i<=x.size()-1;i++){
            res=res+tem*diff_table[i][i];
            tem=tem*Polynomial(-x[i],1);
        }
        return res;
    }
};

#else
#endif
